Load Data

dataset <- read.delim("raw_data/SFigure3_H2O2.txt", stringsAsFactors = FALSE)

dataset$genotype <- gsub(" ", "", dataset$genotype )
dataset$genotype <- factor(dataset$genotype)
dataset$Experiment <- factor(rep(paste0("exp", 1:(length(dataset$genotype)/length(levels(dataset$genotype)))), each=length(unique(dataset$genotype))))

dataset$siRNA <-  factor(gsub(".*[T,O]\\+","",dataset$genotype))
dataset$genotype <-  factor(gsub("\\+.*","",dataset$genotype))

dataset$UID <- factor(paste(dataset$Experiment, dataset$genotype, dataset$siRNA))
dataset$GSID <- factor(paste(dataset$genotype, dataset$siRNA))

# wide format
kable(dataset, row.names = F)
genotype NT H2O2_5uM H2O2_30uM H2O2_50uM Experiment siRNA UID GSID
WT 2420 1528 977 248 exp1 siCtrl exp1 WT siCtrl WT siCtrl
WT 2180 755 303 60 exp1 siPARP1 exp1 WT siPARP1 WT siPARP1
ALC1KO 2428 923 244 95 exp1 siCtrl exp1 ALC1KO siCtrl ALC1KO siCtrl
ALC1KO 2188 834 296 74 exp1 siPARP1 exp1 ALC1KO siPARP1 ALC1KO siPARP1
WT 2072 1332 598 194 exp2 siCtrl exp2 WT siCtrl WT siCtrl
WT 2504 624 344 111 exp2 siPARP1 exp2 WT siPARP1 WT siPARP1
ALC1KO 3160 1222 388 122 exp2 siCtrl exp2 ALC1KO siCtrl ALC1KO siCtrl
ALC1KO 3536 1296 396 33 exp2 siPARP1 exp2 ALC1KO siPARP1 ALC1KO siPARP1
WT 1768 896 562 219 exp3 siCtrl exp3 WT siCtrl WT siCtrl
WT 1888 696 422 160 exp3 siPARP1 exp3 WT siPARP1 WT siPARP1
ALC1KO 2455 712 499 140 exp3 siCtrl exp3 ALC1KO siCtrl ALC1KO siCtrl
ALC1KO 3222 1224 467 67 exp3 siPARP1 exp3 ALC1KO siPARP1 ALC1KO siPARP1
library(reshape2)
# reshape to long format
dataset <- melt(dataset, variable.name = "Treatment", value.name = "Counts")

dataset$genotype <- relevel(dataset$genotype, ref = "WT")
dataset$siRNA <- relevel(dataset$siRNA, ref = "siCtrl")
dataset$UID <- relevel(dataset$UID, ref = "exp1 WT siCtrl")

dataset$H2O2 <- gsub("NT","1",dataset$Treatment)
dataset$H2O2 <- gsub("H2O2_|uM","",dataset$H2O2)
dataset$H2O2 <- log10(as.integer(dataset$H2O2))




dataset$Offset <- NA
for(uid in levels(dataset$UID)){
        dataset$Offset[dataset$UID == uid] <- mean(dataset$Counts[dataset$UID == uid])
}

dataset$NormCounts <- dataset$Counts / dataset$Offset



dataset$Offset2 <- NA
for(gsid in levels(dataset$GSID)){
        dataset$Offset2[dataset$GSID == gsid] <- mean(dataset$NormCounts[dataset$GSID == gsid & dataset$H2O2 == 0])
}

dataset$NormCounts2 <- dataset$NormCounts / dataset$Offset2



# long format
kable(dataset, row.names = F)
genotype Experiment siRNA UID GSID Treatment Counts H2O2 Offset NormCounts Offset2 NormCounts2
WT exp1 siCtrl exp1 WT siCtrl WT siCtrl NT 2420 0.000000 1293.25 1.8712546 1.966433 0.9515984
WT exp1 siPARP1 exp1 WT siPARP1 WT siPARP1 NT 2180 0.000000 824.50 2.6440267 2.608265 1.0137111
ALC1KO exp1 siCtrl exp1 ALC1KO siCtrl ALC1KO siCtrl NT 2428 0.000000 922.50 2.6319783 2.598642 1.0128285
ALC1KO exp1 siPARP1 exp1 ALC1KO siPARP1 ALC1KO siPARP1 NT 2188 0.000000 848.00 2.5801887 2.618868 0.9852307
WT exp2 siCtrl exp2 WT siCtrl WT siCtrl NT 2072 0.000000 1049.00 1.9752145 1.966433 1.0044656
WT exp2 siPARP1 exp2 WT siPARP1 WT siPARP1 NT 2504 0.000000 895.75 2.7954228 2.608265 1.0717558
ALC1KO exp2 siCtrl exp2 ALC1KO siCtrl ALC1KO siCtrl NT 3160 0.000000 1223.00 2.5838103 2.598642 0.9942926
ALC1KO exp2 siPARP1 exp2 ALC1KO siPARP1 ALC1KO siPARP1 NT 3536 0.000000 1315.25 2.6884623 2.618868 1.0265743
WT exp3 siCtrl exp3 WT siCtrl WT siCtrl NT 1768 0.000000 861.25 2.0528302 1.966433 1.0439359
WT exp3 siPARP1 exp3 WT siPARP1 WT siPARP1 NT 1888 0.000000 791.50 2.3853443 2.608265 0.9145331
ALC1KO exp3 siCtrl exp3 ALC1KO siCtrl ALC1KO siCtrl NT 2455 0.000000 951.50 2.5801366 2.598642 0.9928789
ALC1KO exp3 siPARP1 exp3 ALC1KO siPARP1 ALC1KO siPARP1 NT 3222 0.000000 1245.00 2.5879518 2.618868 0.9881950
WT exp1 siCtrl exp1 WT siCtrl WT siCtrl H2O2_5uM 1528 0.698970 1293.25 1.1815194 1.966433 0.6008440
WT exp1 siPARP1 exp1 WT siPARP1 WT siPARP1 H2O2_5uM 755 0.698970 824.50 0.9157065 2.608265 0.3510788
ALC1KO exp1 siCtrl exp1 ALC1KO siCtrl ALC1KO siCtrl H2O2_5uM 923 0.698970 922.50 1.0005420 2.598642 0.3850250
ALC1KO exp1 siPARP1 exp1 ALC1KO siPARP1 ALC1KO siPARP1 H2O2_5uM 834 0.698970 848.00 0.9834906 2.618868 0.3755404
WT exp2 siCtrl exp2 WT siCtrl WT siCtrl H2O2_5uM 1332 0.698970 1049.00 1.2697807 1.966433 0.6457279
WT exp2 siPARP1 exp2 WT siPARP1 WT siPARP1 H2O2_5uM 624 0.698970 895.75 0.6966229 2.608265 0.2670829
ALC1KO exp2 siCtrl exp2 ALC1KO siCtrl ALC1KO siCtrl H2O2_5uM 1222 0.698970 1223.00 0.9991823 2.598642 0.3845018
ALC1KO exp2 siPARP1 exp2 ALC1KO siPARP1 ALC1KO siPARP1 H2O2_5uM 1296 0.698970 1315.25 0.9853640 2.618868 0.3762558
WT exp3 siCtrl exp3 WT siCtrl WT siCtrl H2O2_5uM 896 0.698970 861.25 1.0403483 1.966433 0.5290535
WT exp3 siPARP1 exp3 WT siPARP1 WT siPARP1 H2O2_5uM 696 0.698970 791.50 0.8793430 2.608265 0.3371372
ALC1KO exp3 siCtrl exp3 ALC1KO siCtrl ALC1KO siCtrl H2O2_5uM 712 0.698970 951.50 0.7482922 2.598642 0.2879551
ALC1KO exp3 siPARP1 exp3 ALC1KO siPARP1 ALC1KO siPARP1 H2O2_5uM 1224 0.698970 1245.00 0.9831325 2.618868 0.3754037
WT exp1 siCtrl exp1 WT siCtrl WT siCtrl H2O2_30uM 977 1.477121 1293.25 0.7554610 1.966433 0.3841784
WT exp1 siPARP1 exp1 WT siPARP1 WT siPARP1 H2O2_30uM 303 1.477121 824.50 0.3674955 2.608265 0.1408965
ALC1KO exp1 siCtrl exp1 ALC1KO siCtrl ALC1KO siCtrl H2O2_30uM 244 1.477121 922.50 0.2644986 2.598642 0.1017834
ALC1KO exp1 siPARP1 exp1 ALC1KO siPARP1 ALC1KO siPARP1 H2O2_30uM 296 1.477121 848.00 0.3490566 2.618868 0.1332853
WT exp2 siCtrl exp2 WT siCtrl WT siCtrl H2O2_30uM 598 1.477121 1049.00 0.5700667 1.966433 0.2898989
WT exp2 siPARP1 exp2 WT siPARP1 WT siPARP1 H2O2_30uM 344 1.477121 895.75 0.3840357 2.608265 0.1472380
ALC1KO exp2 siCtrl exp2 ALC1KO siCtrl ALC1KO siCtrl H2O2_30uM 388 1.477121 1223.00 0.3172527 2.598642 0.1220840
ALC1KO exp2 siPARP1 exp2 ALC1KO siPARP1 ALC1KO siPARP1 H2O2_30uM 396 1.477121 1315.25 0.3010834 2.618868 0.1149670
WT exp3 siCtrl exp3 WT siCtrl WT siCtrl H2O2_30uM 562 1.477121 861.25 0.6525399 1.966433 0.3318394
WT exp3 siPARP1 exp3 WT siPARP1 WT siPARP1 H2O2_30uM 422 1.477121 791.50 0.5331649 2.608265 0.2044136
ALC1KO exp3 siCtrl exp3 ALC1KO siCtrl ALC1KO siCtrl H2O2_30uM 499 1.477121 951.50 0.5244351 2.598642 0.2018112
ALC1KO exp3 siPARP1 exp3 ALC1KO siPARP1 ALC1KO siPARP1 H2O2_30uM 467 1.477121 1245.00 0.3751004 2.618868 0.1432300
WT exp1 siCtrl exp1 WT siCtrl WT siCtrl H2O2_50uM 248 1.698970 1293.25 0.1917649 1.966433 0.0975192
WT exp1 siPARP1 exp1 WT siPARP1 WT siPARP1 H2O2_50uM 60 1.698970 824.50 0.0727714 2.608265 0.0279003
ALC1KO exp1 siCtrl exp1 ALC1KO siCtrl ALC1KO siCtrl H2O2_50uM 95 1.698970 922.50 0.1029810 2.598642 0.0396288
ALC1KO exp1 siPARP1 exp1 ALC1KO siPARP1 ALC1KO siPARP1 H2O2_50uM 74 1.698970 848.00 0.0872642 2.618868 0.0333213
WT exp2 siCtrl exp2 WT siCtrl WT siCtrl H2O2_50uM 194 1.698970 1049.00 0.1849380 1.966433 0.0940475
WT exp2 siPARP1 exp2 WT siPARP1 WT siPARP1 H2O2_50uM 111 1.698970 895.75 0.1239185 2.608265 0.0475099
ALC1KO exp2 siCtrl exp2 ALC1KO siCtrl ALC1KO siCtrl H2O2_50uM 122 1.698970 1223.00 0.0997547 2.598642 0.0383872
ALC1KO exp2 siPARP1 exp2 ALC1KO siPARP1 ALC1KO siPARP1 H2O2_50uM 33 1.698970 1315.25 0.0250903 2.618868 0.0095806
WT exp3 siCtrl exp3 WT siCtrl WT siCtrl H2O2_50uM 219 1.698970 861.25 0.2542816 1.966433 0.1293111
WT exp3 siPARP1 exp3 WT siPARP1 WT siPARP1 H2O2_50uM 160 1.698970 791.50 0.2021478 2.608265 0.0775028
ALC1KO exp3 siCtrl exp3 ALC1KO siCtrl ALC1KO siCtrl H2O2_50uM 140 1.698970 951.50 0.1471361 2.598642 0.0566204
ALC1KO exp3 siPARP1 exp3 ALC1KO siPARP1 ALC1KO siPARP1 H2O2_50uM 67 1.698970 1245.00 0.0538153 2.618868 0.0205491

Plot Data

library(ggplot2)

# raw data
ggplot(dataset, aes(x=H2O2, y=Counts)) + 
        theme_bw() +
        theme(panel.grid=element_blank(), text = element_text(size=14)) +
        geom_smooth(method=lm, formula = y ~ poly(x,2), se=FALSE, aes(colour=siRNA)) +
        geom_point(aes(colour=siRNA, shape=Experiment), size=2) +        
        facet_grid(. ~ genotype) +
        xlab(label = "H2O2 (log10 µM)") +
        scale_shape_manual(values=15:20) +
        scale_color_manual(values=c("#000000","#FF0000"))

# NormCounts Linear
ggplot(dataset, aes(x=H2O2, y=NormCounts, color=siRNA)) + 
        theme_bw() +
        theme(panel.grid=element_blank(), text = element_text(size=14)) +
        geom_point(aes(colour=siRNA), size=2) +        
        geom_smooth(method=lm, formula = y ~ x, se=FALSE) +
        facet_grid(. ~ genotype) +
        xlab(label = "H2O2 (log10 µM)") +
        scale_color_manual(values=c("#000000","#FF0000"))

# NormCounts2 Linear
ggplot(dataset, aes(x=H2O2, y=NormCounts2, color=siRNA)) + 
        theme_bw() +
        theme(panel.grid=element_blank(), text = element_text(size=14)) +
        geom_point(aes(colour=siRNA), size=2) +        
        geom_smooth(method=lm, formula = y ~ x, se=FALSE) +
        facet_grid(. ~ genotype) +
        xlab(label = "H2O2 (log10 µM)") +
        scale_color_manual(values=c("#000000","#FF0000"))

# NormCounts Quadratic
ggplot(dataset, aes(x=H2O2, y=NormCounts, color=siRNA)) + 
        theme_bw() +
        theme(panel.grid=element_blank(), text = element_text(size=14)) +
        geom_point(aes(colour=siRNA), size=2) +        
        geom_smooth(method=lm, formula = y ~ poly(x,2), se=FALSE) +
        facet_grid(. ~ genotype) +
        xlab(label = "H2O2 (log10 µM)")+
        scale_color_manual(values=c("#000000","#FF0000"))

# NormCounts2 Quadratic
ggplot(dataset, aes(x=H2O2, y=NormCounts2, color=siRNA)) + 
        theme_bw() +
        theme(panel.grid=element_blank(), text = element_text(size=14)) +
        geom_point(aes(colour=siRNA), size=2) +        
        geom_smooth(method=lm, formula = y ~ poly(x,2), se=FALSE) +
        facet_grid(. ~ genotype) +
        xlab(label = "H2O2 (log10 µM)") +
        scale_color_manual(values=c("#000000","#FF0000"))

# NormCounts Cubic
ggplot(dataset, aes(x=H2O2, y=NormCounts, color=siRNA)) + 
        theme_bw() +
        theme(panel.grid=element_blank(), text = element_text(size=14)) +
        geom_point(aes(colour=siRNA), size=2) +        
        geom_smooth(method=lm, formula = y ~ poly(x,3), se=FALSE) +
        facet_grid(. ~ genotype) +
        xlab(label = "H2O2 (log10 µM)")+
        scale_color_manual(values=c("#000000","#FF0000"))

# NormCounts2 Cubic
ggplot(dataset, aes(x=H2O2, y=NormCounts2, color=siRNA)) + 
        theme_bw() +
        theme(panel.grid=element_blank(), text = element_text(size=14)) +
        geom_point(aes(colour=siRNA), size=2) +        
        geom_smooth(method=lm, formula = y ~ poly(x,3), se=FALSE) +
        facet_grid(. ~ genotype) +
        xlab(label = "H2O2 (log10 µM)") +
        scale_color_manual(values=c("#000000","#FF0000"))

pdf("FigureS3_H2O2.pdf", width = 5, height = 4)

ggplot(dataset, aes(x=H2O2, y=NormCounts2)) + 
        theme_bw() +
        theme(panel.grid=element_blank(), text = element_text(size=14)) +
        geom_point(aes(colour = siRNA, shape = genotype), size=1.75) +
        geom_smooth(method=lm, formula = y ~ poly(x,2), se=TRUE, 
                    aes(group = GSID,colour = siRNA, linetype = genotype), fill='#DDDDDD', size=0.5) +
        #facet_grid(. ~ genotype) +
        xlab(label = "H2O2 (log10 µM)") +
        ylab(label = "Normalized Counts") +
        scale_color_manual(values=c("#000000","#FF0000")) +
        guides(linetype = guide_legend(override.aes= list(color = "#555555"))) 

dev.off()
## quartz_off_screen 
##                 2

Models

library(MASS)
library(DHARMa)
library(lme4)
library(lmerTest)
library(bbmle)

Linear formula

fit1 <- lm(Counts ~ Experiment + H2O2*siRNA*genotype, data = dataset)
print(summary(fit1))
## 
## Call:
## lm(formula = Counts ~ Experiment + H2O2 * siRNA * genotype, data = dataset)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -697.48 -230.97    3.66  207.16  711.48 
## 
## Coefficients:
##                                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                       2008.20     197.38  10.174 2.11e-12 ***
## Experimentexp2                     148.69     127.73   1.164   0.2516    
## Experimentexp3                      -9.75     127.73  -0.076   0.9396    
## H2O2                             -1018.50     155.34  -6.557 9.85e-08 ***
## siRNAsiPARP1                      -117.27     258.93  -0.453   0.6532    
## genotypeALC1KO                     361.79     258.93   1.397   0.1704    
## H2O2:siRNAsiPARP1                 -116.96     219.68  -0.532   0.5975    
## H2O2:genotypeALC1KO               -410.10     219.68  -1.867   0.0697 .  
## siRNAsiPARP1:genotypeALC1KO        423.11     366.19   1.155   0.2551    
## H2O2:siRNAsiPARP1:genotypeALC1KO   -91.64     310.68  -0.295   0.7696    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 361.3 on 38 degrees of freedom
## Multiple R-squared:  0.8871, Adjusted R-squared:  0.8603 
## F-statistic: 33.17 on 9 and 38 DF,  p-value: 2.56e-15
cat("AIC: ", AIC(fit1))
## AIC:  712.4085
simres <- simulateResiduals(fittedModel = fit1)
plot(simres)

fit2 <- lm(NormCounts ~ H2O2*siRNA*genotype, data = dataset)
print(summary(fit2))
## 
## Call:
## lm(formula = NormCounts ~ H2O2 * siRNA * genotype, data = dataset)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.66792 -0.08492  0.07217  0.16452  0.48643 
## 
## Coefficients:
##                                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                        1.9297     0.1435  13.447  < 2e-16 ***
## H2O2                              -0.9597     0.1218  -7.882 1.14e-09 ***
## siRNAsiPARP1                       0.3793     0.2029   1.869   0.0690 .  
## genotypeALC1KO                     0.4093     0.2029   2.017   0.0505 .  
## H2O2:siRNAsiPARP1                 -0.3915     0.1722  -2.274   0.0284 *  
## H2O2:genotypeALC1KO               -0.4225     0.1722  -2.454   0.0186 *  
## siRNAsiPARP1:genotypeALC1KO       -0.3275     0.2870  -1.141   0.2606    
## H2O2:siRNAsiPARP1:genotypeALC1KO   0.3381     0.2435   1.388   0.1727    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2832 on 40 degrees of freedom
## Multiple R-squared:  0.9189, Adjusted R-squared:  0.9047 
## F-statistic: 64.75 on 7 and 40 DF,  p-value: < 2.2e-16
cat("AIC: ", AIC(fit2))
## AIC:  24.33749
simres <- simulateResiduals(fittedModel = fit2)
plot(simres)

fit3 <- lm(NormCounts2 ~ H2O2*siRNA*genotype, data = dataset)
print(summary(fit3))
## 
## Call:
## lm(formula = NormCounts2 ~ H2O2 * siRNA * genotype, data = dataset)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.25608 -0.04318  0.02865  0.07163  0.18649 
## 
## Coefficients:
##                                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                       0.98133    0.05621  17.457  < 2e-16 ***
## H2O2                             -0.48804    0.04769 -10.233  9.9e-13 ***
## siRNAsiPARP1                     -0.09607    0.07950  -1.208    0.234    
## genotypeALC1KO                   -0.08123    0.07950  -1.022    0.313    
## H2O2:siRNAsiPARP1                -0.03001    0.06745  -0.445    0.659    
## H2O2:genotypeALC1KO              -0.04386    0.06745  -0.650    0.519    
## siRNAsiPARP1:genotypeALC1KO       0.10889    0.11243   0.969    0.339    
## H2O2:siRNAsiPARP1:genotypeALC1KO  0.01371    0.09539   0.144    0.886    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1109 on 40 degrees of freedom
## Multiple R-squared:  0.9246, Adjusted R-squared:  0.9114 
## F-statistic: 70.09 on 7 and 40 DF,  p-value: < 2.2e-16
cat("AIC: ", AIC(fit3))
## AIC:  -65.63346
simres <- simulateResiduals(fittedModel = fit3)
plot(simres)

fit4 <- lmer(Counts ~ H2O2*siRNA*genotype + (1|UID), data = dataset)
print(summary(fit4))
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Counts ~ H2O2 * siRNA * genotype + (1 | UID)
##    Data: dataset
## 
## REML criterion at convergence: 601.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.0069 -0.4188  0.1341  0.4384  2.2311 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  UID      (Intercept)   3255    57.06  
##  Residual             127693   357.34  
## Number of obs: 48, groups:  UID, 12
## 
## Fixed effects:
##                                  Estimate Std. Error       df t value Pr(>|t|)
## (Intercept)                       2054.52     184.07    35.29  11.161 3.97e-13
## H2O2                             -1018.50     153.65    32.00  -6.629 1.78e-07
## siRNAsiPARP1                      -117.27     260.32    35.29  -0.451   0.6551
## genotypeALC1KO                     361.79     260.32    35.29   1.390   0.1733
## H2O2:siRNAsiPARP1                 -116.96     217.29    32.00  -0.538   0.5941
## H2O2:genotypeALC1KO               -410.10     217.29    32.00  -1.887   0.0682
## siRNAsiPARP1:genotypeALC1KO        423.11     368.15    35.29   1.149   0.2582
## H2O2:siRNAsiPARP1:genotypeALC1KO   -91.64     307.30    32.00  -0.298   0.7675
##                                     
## (Intercept)                      ***
## H2O2                             ***
## siRNAsiPARP1                        
## genotypeALC1KO                      
## H2O2:siRNAsiPARP1                   
## H2O2:genotypeALC1KO              .  
## siRNAsiPARP1:genotypeALC1KO         
## H2O2:siRNAsiPARP1:genotypeALC1KO    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##                (Intr) H2O2   sRNAsPARP1 gALC1K H2O2:sRNAPARP1 H2O2:A sRNAPARP1:
## H2O2           -0.809                                                          
## siRNAsPARP1    -0.707  0.572                                                   
## gntypALC1KO    -0.707  0.572  0.500                                            
## H2O2:sRNAPARP1  0.572 -0.707 -0.809     -0.404                                 
## H2O2:ALC1KO     0.572 -0.707 -0.404     -0.809  0.500                          
## sRNAPARP1:A     0.500 -0.404 -0.707     -0.707  0.572          0.572           
## H2O2:RNAPARP1: -0.404  0.500  0.572      0.572 -0.707         -0.707 -0.809
cat("AIC: ", AIC(fit4))
## AIC:  621.279
simres <- simulateResiduals(fittedModel = fit4)
plot(simres)

Quadratic formula

fit5 <- lm(Counts ~ Experiment + poly(H2O2,2)*siRNA*genotype, data = dataset)
print(summary(fit5))
## 
## Call:
## lm(formula = Counts ~ Experiment + poly(H2O2, 2) * siRNA * genotype, 
##     data = dataset)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -725.74 -158.61  -17.09  102.58  475.60 
## 
## Coefficients:
##                                            Estimate Std. Error t value Pr(>|t|)
## (Intercept)                                 1021.52      98.18  10.405 4.18e-12
## Experimentexp2                               148.69      98.18   1.514   0.1392
## Experimentexp3                                -9.75      98.18  -0.099   0.9215
## poly(H2O2, 2)1                             -4737.38     555.39  -8.530 5.80e-10
## poly(H2O2, 2)2                                61.12     555.39   0.110   0.9130
## siRNAsiPARP1                                -230.58     113.37  -2.034   0.0498
## genotypeALC1KO                               -35.50     113.37  -0.313   0.7561
## poly(H2O2, 2)1:siRNAsiPARP1                 -544.03     785.44  -0.693   0.4932
## poly(H2O2, 2)2:siRNAsiPARP1                 1647.83     785.44   2.098   0.0434
## poly(H2O2, 2)1:genotypeALC1KO              -1907.53     785.44  -2.429   0.0206
## poly(H2O2, 2)2:genotypeALC1KO               1754.59     785.44   2.234   0.0322
## siRNAsiPARP1:genotypeALC1KO                  334.33     160.33   2.085   0.0446
## poly(H2O2, 2)1:siRNAsiPARP1:genotypeALC1KO  -426.25    1110.78  -0.384   0.7036
## poly(H2O2, 2)2:siRNAsiPARP1:genotypeALC1KO -1694.41    1110.78  -1.525   0.1364
##                                               
## (Intercept)                                ***
## Experimentexp2                                
## Experimentexp3                                
## poly(H2O2, 2)1                             ***
## poly(H2O2, 2)2                                
## siRNAsiPARP1                               *  
## genotypeALC1KO                                
## poly(H2O2, 2)1:siRNAsiPARP1                   
## poly(H2O2, 2)2:siRNAsiPARP1                *  
## poly(H2O2, 2)1:genotypeALC1KO              *  
## poly(H2O2, 2)2:genotypeALC1KO              *  
## siRNAsiPARP1:genotypeALC1KO                *  
## poly(H2O2, 2)1:siRNAsiPARP1:genotypeALC1KO    
## poly(H2O2, 2)2:siRNAsiPARP1:genotypeALC1KO    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 277.7 on 34 degrees of freedom
## Multiple R-squared:  0.9403, Adjusted R-squared:  0.9175 
## F-statistic:  41.2 on 13 and 34 DF,  p-value: < 2.2e-16
cat("AIC: ", AIC(fit5))
## AIC:  689.8121
simres <- simulateResiduals(fittedModel = fit5)
plot(simres)

fit6 <- lm(NormCounts ~ poly(H2O2,2)*siRNA*genotype, data = dataset)
print(summary(fit6))
## 
## Call:
## lm(formula = NormCounts ~ poly(H2O2, 2) * siRNA * genotype, data = dataset)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.233442 -0.077535 -0.008586  0.088724  0.298066 
## 
## Coefficients:
##                                              Estimate Std. Error t value
## (Intercept)                                 1.000e+00  4.284e-02  23.343
## poly(H2O2, 2)1                             -4.464e+00  2.968e-01 -15.040
## poly(H2O2, 2)2                              1.093e-01  2.968e-01   0.368
## siRNAsiPARP1                                7.412e-17  6.058e-02   0.000
## genotypeALC1KO                             -5.228e-17  6.058e-02   0.000
## poly(H2O2, 2)1:siRNAsiPARP1                -1.821e+00  4.197e-01  -4.338
## poly(H2O2, 2)2:siRNAsiPARP1                 1.907e+00  4.197e-01   4.544
## poly(H2O2, 2)1:genotypeALC1KO              -1.965e+00  4.197e-01  -4.682
## poly(H2O2, 2)2:genotypeALC1KO               1.670e+00  4.197e-01   3.978
## siRNAsiPARP1:genotypeALC1KO                 9.821e-17  8.568e-02   0.000
## poly(H2O2, 2)1:siRNAsiPARP1:genotypeALC1KO  1.572e+00  5.936e-01   2.649
## poly(H2O2, 2)2:siRNAsiPARP1:genotypeALC1KO -2.133e+00  5.936e-01  -3.593
##                                            Pr(>|t|)    
## (Intercept)                                 < 2e-16 ***
## poly(H2O2, 2)1                              < 2e-16 ***
## poly(H2O2, 2)2                             0.714841    
## siRNAsiPARP1                               1.000000    
## genotypeALC1KO                             1.000000    
## poly(H2O2, 2)1:siRNAsiPARP1                0.000111 ***
## poly(H2O2, 2)2:siRNAsiPARP1                6.00e-05 ***
## poly(H2O2, 2)1:genotypeALC1KO              3.95e-05 ***
## poly(H2O2, 2)2:genotypeALC1KO              0.000322 ***
## siRNAsiPARP1:genotypeALC1KO                1.000000    
## poly(H2O2, 2)1:siRNAsiPARP1:genotypeALC1KO 0.011915 *  
## poly(H2O2, 2)2:siRNAsiPARP1:genotypeALC1KO 0.000970 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1484 on 36 degrees of freedom
## Multiple R-squared:   0.98,  Adjusted R-squared:  0.9738 
## F-statistic:   160 on 11 and 36 DF,  p-value: < 2.2e-16
cat("AIC: ", AIC(fit6))
## AIC:  -34.74272
simres <- simulateResiduals(fittedModel = fit6)
plot(simres)

fit7 <- lm(NormCounts2 ~ poly(H2O2,2)*siRNA*genotype, data = dataset)
print(summary(fit7))
## 
## Call:
## lm(formula = NormCounts2 ~ poly(H2O2, 2) * siRNA * genotype, 
##     data = dataset)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.098930 -0.035263 -0.003283  0.033948  0.125891 
## 
## Coefficients:
##                                            Estimate Std. Error t value Pr(>|t|)
## (Intercept)                                 0.50853    0.01779  28.581  < 2e-16
## poly(H2O2, 2)1                             -2.27004    0.12327 -18.415  < 2e-16
## poly(H2O2, 2)2                              0.05558    0.12327   0.451 0.654770
## siRNAsiPARP1                               -0.12514    0.02516  -4.973 1.63e-05
## genotypeALC1KO                             -0.12372    0.02516  -4.917 1.94e-05
## poly(H2O2, 2)1:siRNAsiPARP1                -0.13957    0.17433  -0.801 0.428600
## poly(H2O2, 2)2:siRNAsiPARP1                 0.71759    0.17433   4.116 0.000214
## poly(H2O2, 2)1:genotypeALC1KO              -0.20399    0.17433  -1.170 0.249633
## poly(H2O2, 2)2:genotypeALC1KO               0.62899    0.17433   3.608 0.000930
## siRNAsiPARP1:genotypeALC1KO                 0.12217    0.03559   3.433 0.001517
## poly(H2O2, 2)1:siRNAsiPARP1:genotypeALC1KO  0.06375    0.24654   0.259 0.797440
## poly(H2O2, 2)2:siRNAsiPARP1:genotypeALC1KO -0.80896    0.24654  -3.281 0.002301
##                                               
## (Intercept)                                ***
## poly(H2O2, 2)1                             ***
## poly(H2O2, 2)2                                
## siRNAsiPARP1                               ***
## genotypeALC1KO                             ***
## poly(H2O2, 2)1:siRNAsiPARP1                   
## poly(H2O2, 2)2:siRNAsiPARP1                ***
## poly(H2O2, 2)1:genotypeALC1KO                 
## poly(H2O2, 2)2:genotypeALC1KO              ***
## siRNAsiPARP1:genotypeALC1KO                ** 
## poly(H2O2, 2)1:siRNAsiPARP1:genotypeALC1KO    
## poly(H2O2, 2)2:siRNAsiPARP1:genotypeALC1KO ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.06164 on 36 degrees of freedom
## Multiple R-squared:  0.9791, Adjusted R-squared:  0.9727 
## F-statistic:   153 on 11 and 36 DF,  p-value: < 2.2e-16
cat("AIC: ", AIC(fit7))
## AIC:  -119.0964
simres <- simulateResiduals(fittedModel = fit7)
plot(simres)

fit8 <- lmer(Counts ~ poly(H2O2,2)*siRNA*genotype + (1|UID), data = dataset)
print(summary(fit8))
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Counts ~ poly(H2O2, 2) * siRNA * genotype + (1 | UID)
##    Data: dataset
## 
## REML criterion at convergence: 505
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.44832 -0.62194  0.01954  0.46624  1.90599 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  UID      (Intercept) 19567    139.9   
##  Residual             62446    249.9   
## Number of obs: 48, groups:  UID, 12
## 
## Fixed effects:
##                                            Estimate Std. Error       df t value
## (Intercept)                                 1067.83     108.29     8.00   9.861
## poly(H2O2, 2)1                             -4737.38     499.78    28.00  -9.479
## poly(H2O2, 2)2                                61.12     499.78    28.00   0.122
## siRNAsiPARP1                                -230.58     153.14     8.00  -1.506
## genotypeALC1KO                               -35.50     153.14     8.00  -0.232
## poly(H2O2, 2)1:siRNAsiPARP1                 -544.03     706.80    28.00  -0.770
## poly(H2O2, 2)2:siRNAsiPARP1                 1647.83     706.80    28.00   2.331
## poly(H2O2, 2)1:genotypeALC1KO              -1907.53     706.80    28.00  -2.699
## poly(H2O2, 2)2:genotypeALC1KO               1754.59     706.80    28.00   2.482
## siRNAsiPARP1:genotypeALC1KO                  334.33     216.57     8.00   1.544
## poly(H2O2, 2)1:siRNAsiPARP1:genotypeALC1KO  -426.25     999.57    28.00  -0.426
## poly(H2O2, 2)2:siRNAsiPARP1:genotypeALC1KO -1694.41     999.57    28.00  -1.695
##                                            Pr(>|t|)    
## (Intercept)                                9.42e-06 ***
## poly(H2O2, 2)1                             3.09e-10 ***
## poly(H2O2, 2)2                               0.9035    
## siRNAsiPARP1                                 0.1706    
## genotypeALC1KO                               0.8225    
## poly(H2O2, 2)1:siRNAsiPARP1                  0.4479    
## poly(H2O2, 2)2:siRNAsiPARP1                  0.0272 *  
## poly(H2O2, 2)1:genotypeALC1KO                0.0117 *  
## poly(H2O2, 2)2:genotypeALC1KO                0.0193 *  
## siRNAsiPARP1:genotypeALC1KO                  0.1612    
## poly(H2O2, 2)1:siRNAsiPARP1:genotypeALC1KO   0.6731    
## poly(H2O2, 2)2:siRNAsiPARP1:genotypeALC1KO   0.1011    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##                      (Intr) pl(H2O2,2)1 pl(H2O2,2)2 sRNAsPARP1 gALC1K
## pl(H2O2,2)1           0.000                                          
## pl(H2O2,2)2           0.000  0.000                                   
## siRNAsPARP1          -0.707  0.000       0.000                       
## gntypALC1KO          -0.707  0.000       0.000       0.500           
## pl(H2O2,2)1:RNAPARP1  0.000 -0.707       0.000       0.000      0.000
## pl(H2O2,2)2:RNAPARP1  0.000  0.000      -0.707       0.000      0.000
## p(H2O2,2)1:A          0.000 -0.707       0.000       0.000      0.000
## p(H2O2,2)2:A          0.000  0.000      -0.707       0.000      0.000
## sRNAPARP1:A           0.500  0.000       0.000      -0.707     -0.707
## p(H2O2,2)1:RNAPARP1:  0.000  0.500       0.000       0.000      0.000
## p(H2O2,2)2:RNAPARP1:  0.000  0.000       0.500       0.000      0.000
##                      pl(H2O2,2)1:RNAPARP1 pl(H2O2,2)2:RNAPARP1 p(H2O2,2)1:A
## pl(H2O2,2)1                                                                
## pl(H2O2,2)2                                                                
## siRNAsPARP1                                                                
## gntypALC1KO                                                                
## pl(H2O2,2)1:RNAPARP1                                                       
## pl(H2O2,2)2:RNAPARP1  0.000                                                
## p(H2O2,2)1:A          0.500                0.000                           
## p(H2O2,2)2:A          0.000                0.500                0.000      
## sRNAPARP1:A           0.000                0.000                0.000      
## p(H2O2,2)1:RNAPARP1: -0.707                0.000               -0.707      
## p(H2O2,2)2:RNAPARP1:  0.000               -0.707                0.000      
##                      p(H2O2,2)2:A sRNAPARP1: p(H2O2,2)1:RNAPARP1:
## pl(H2O2,2)1                                                      
## pl(H2O2,2)2                                                      
## siRNAsPARP1                                                      
## gntypALC1KO                                                      
## pl(H2O2,2)1:RNAPARP1                                             
## pl(H2O2,2)2:RNAPARP1                                             
## p(H2O2,2)1:A                                                     
## p(H2O2,2)2:A                                                     
## sRNAPARP1:A           0.000                                      
## p(H2O2,2)1:RNAPARP1:  0.000        0.000                         
## p(H2O2,2)2:RNAPARP1: -0.707        0.000      0.000
cat("AIC: ", AIC(fit8))
## AIC:  533.0262
simres <- simulateResiduals(fittedModel = fit8)
plot(simres)

Cubic formula

fit9 <- lm(Counts ~ Experiment + poly(H2O2,3)*siRNA*genotype, data = dataset)
print(summary(fit9))
## 
## Call:
## lm(formula = Counts ~ Experiment + poly(H2O2, 3) * siRNA * genotype, 
##     data = dataset)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -747.69 -119.62    6.48  112.67  451.62 
## 
## Coefficients:
##                                             Estimate Std. Error t value
## (Intercept)                                 1021.521     94.245  10.839
## Experimentexp2                               148.688     94.245   1.578
## Experimentexp3                                -9.750     94.245  -0.103
## poly(H2O2, 3)1                             -4737.379    533.132  -8.886
## poly(H2O2, 3)2                                61.116    533.132   0.115
## poly(H2O2, 3)3                              -743.243    533.132  -1.394
## siRNAsiPARP1                                -230.583    108.825  -2.119
## genotypeALC1KO                               -35.500    108.825  -0.326
## poly(H2O2, 3)1:siRNAsiPARP1                 -544.032    753.963  -0.722
## poly(H2O2, 3)2:siRNAsiPARP1                 1647.828    753.963   2.186
## poly(H2O2, 3)3:siRNAsiPARP1                    9.997    753.963   0.013
## poly(H2O2, 3)1:genotypeALC1KO              -1907.534    753.963  -2.530
## poly(H2O2, 3)2:genotypeALC1KO               1754.587    753.963   2.327
## poly(H2O2, 3)3:genotypeALC1KO                106.013    753.963   0.141
## siRNAsiPARP1:genotypeALC1KO                  334.333    153.902   2.172
## poly(H2O2, 3)1:siRNAsiPARP1:genotypeALC1KO  -426.247   1066.265  -0.400
## poly(H2O2, 3)2:siRNAsiPARP1:genotypeALC1KO -1694.412   1066.265  -1.589
## poly(H2O2, 3)3:siRNAsiPARP1:genotypeALC1KO   -54.356   1066.265  -0.051
##                                            Pr(>|t|)    
## (Intercept)                                6.78e-12 ***
## Experimentexp2                               0.1251    
## Experimentexp3                               0.9183    
## poly(H2O2, 3)1                             6.64e-10 ***
## poly(H2O2, 3)2                               0.9095    
## poly(H2O2, 3)3                               0.1735    
## siRNAsiPARP1                                 0.0425 *  
## genotypeALC1KO                               0.7465    
## poly(H2O2, 3)1:siRNAsiPARP1                  0.4761    
## poly(H2O2, 3)2:siRNAsiPARP1                  0.0368 *  
## poly(H2O2, 3)3:siRNAsiPARP1                  0.9895    
## poly(H2O2, 3)1:genotypeALC1KO                0.0169 *  
## poly(H2O2, 3)2:genotypeALC1KO                0.0269 *  
## poly(H2O2, 3)3:genotypeALC1KO                0.8891    
## siRNAsiPARP1:genotypeALC1KO                  0.0379 *  
## poly(H2O2, 3)1:siRNAsiPARP1:genotypeALC1KO   0.6922    
## poly(H2O2, 3)2:siRNAsiPARP1:genotypeALC1KO   0.1225    
## poly(H2O2, 3)3:siRNAsiPARP1:genotypeALC1KO   0.9597    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 266.6 on 30 degrees of freedom
## Multiple R-squared:  0.9515, Adjusted R-squared:  0.924 
## F-statistic:  34.6 on 17 and 30 DF,  p-value: 4.022e-15
cat("AIC: ", AIC(fit9))
## AIC:  687.8777
simres <- simulateResiduals(fittedModel = fit9)
plot(simres)

fit10 <- lm(NormCounts ~ poly(H2O2,3)*siRNA*genotype, data = dataset)
print(summary(fit10))
## 
## Call:
## lm(formula = NormCounts ~ poly(H2O2, 3) * siRNA * genotype, data = dataset)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.222920 -0.039175 -0.001219  0.045161  0.187158 
## 
## Coefficients:
##                                              Estimate Std. Error t value
## (Intercept)                                 1.000e+00  2.783e-02  35.936
## poly(H2O2, 3)1                             -4.464e+00  1.928e-01 -23.154
## poly(H2O2, 3)2                              1.093e-01  1.928e-01   0.567
## poly(H2O2, 3)3                             -6.838e-01  1.928e-01  -3.547
## siRNAsiPARP1                                1.666e-16  3.935e-02   0.000
## genotypeALC1KO                              1.302e-16  3.935e-02   0.000
## poly(H2O2, 3)1:siRNAsiPARP1                -1.821e+00  2.727e-01  -6.679
## poly(H2O2, 3)2:siRNAsiPARP1                 1.907e+00  2.727e-01   6.996
## poly(H2O2, 3)3:siRNAsiPARP1                -1.882e-01  2.727e-01  -0.690
## poly(H2O2, 3)1:genotypeALC1KO              -1.965e+00  2.727e-01  -7.208
## poly(H2O2, 3)2:genotypeALC1KO               1.670e+00  2.727e-01   6.124
## poly(H2O2, 3)3:genotypeALC1KO               5.081e-02  2.727e-01   0.186
## siRNAsiPARP1:genotypeALC1KO                -3.274e-16  5.565e-02   0.000
## poly(H2O2, 3)1:siRNAsiPARP1:genotypeALC1KO  1.572e+00  3.856e-01   4.078
## poly(H2O2, 3)2:siRNAsiPARP1:genotypeALC1KO -2.133e+00  3.856e-01  -5.531
## poly(H2O2, 3)3:siRNAsiPARP1:genotypeALC1KO  2.267e-01  3.856e-01   0.588
##                                            Pr(>|t|)    
## (Intercept)                                 < 2e-16 ***
## poly(H2O2, 3)1                              < 2e-16 ***
## poly(H2O2, 3)2                             0.574725    
## poly(H2O2, 3)3                             0.001227 ** 
## siRNAsiPARP1                               1.000000    
## genotypeALC1KO                             1.000000    
## poly(H2O2, 3)1:siRNAsiPARP1                1.54e-07 ***
## poly(H2O2, 3)2:siRNAsiPARP1                6.30e-08 ***
## poly(H2O2, 3)3:siRNAsiPARP1                0.495028    
## poly(H2O2, 3)1:genotypeALC1KO              3.47e-08 ***
## poly(H2O2, 3)2:genotypeALC1KO              7.59e-07 ***
## poly(H2O2, 3)3:genotypeALC1KO              0.853350    
## siRNAsiPARP1:genotypeALC1KO                1.000000    
## poly(H2O2, 3)1:siRNAsiPARP1:genotypeALC1KO 0.000281 ***
## poly(H2O2, 3)2:siRNAsiPARP1:genotypeALC1KO 4.23e-06 ***
## poly(H2O2, 3)3:siRNAsiPARP1:genotypeALC1KO 0.560695    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0964 on 32 degrees of freedom
## Multiple R-squared:  0.9925, Adjusted R-squared:  0.989 
## F-statistic: 281.6 on 15 and 32 DF,  p-value: < 2.2e-16
cat("AIC: ", AIC(fit10))
## AIC:  -73.81496
simres <- simulateResiduals(fittedModel = fit10)
plot(simres)

fit11 <- lm(NormCounts2 ~ poly(H2O2,3)*siRNA*genotype, data = dataset)
print(summary(fit11))
## 
## Call:
## lm(formula = NormCounts2 ~ poly(H2O2, 3) * siRNA * genotype, 
##     data = dataset)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.085467 -0.014959 -0.000465  0.019616  0.071756 
## 
## Coefficients:
##                                             Estimate Std. Error t value
## (Intercept)                                 0.508535   0.011511  44.179
## poly(H2O2, 3)1                             -2.270035   0.079750 -28.464
## poly(H2O2, 3)2                              0.055582   0.079750   0.697
## poly(H2O2, 3)3                             -0.347722   0.079750  -4.360
## siRNAsiPARP1                               -0.125138   0.016279  -7.687
## genotypeALC1KO                             -0.123719   0.016279  -7.600
## poly(H2O2, 3)1:siRNAsiPARP1                -0.139574   0.112783  -1.238
## poly(H2O2, 3)2:siRNAsiPARP1                 0.717592   0.112783   6.363
## poly(H2O2, 3)3:siRNAsiPARP1                 0.013414   0.112783   0.119
## poly(H2O2, 3)1:genotypeALC1KO              -0.203992   0.112783  -1.809
## poly(H2O2, 3)2:genotypeALC1KO               0.628990   0.112783   5.577
## poly(H2O2, 3)3:genotypeALC1KO               0.104147   0.112783   0.923
## siRNAsiPARP1:genotypeALC1KO                 0.122166   0.023022   5.307
## poly(H2O2, 3)1:siRNAsiPARP1:genotypeALC1KO  0.063748   0.159500   0.400
## poly(H2O2, 3)2:siRNAsiPARP1:genotypeALC1KO -0.808958   0.159500  -5.072
## poly(H2O2, 3)3:siRNAsiPARP1:genotypeALC1KO  0.003173   0.159500   0.020
##                                            Pr(>|t|)    
## (Intercept)                                 < 2e-16 ***
## poly(H2O2, 3)1                              < 2e-16 ***
## poly(H2O2, 3)2                             0.490862    
## poly(H2O2, 3)3                             0.000126 ***
## siRNAsiPARP1                               9.20e-09 ***
## genotypeALC1KO                             1.17e-08 ***
## poly(H2O2, 3)1:siRNAsiPARP1                0.224889    
## poly(H2O2, 3)2:siRNAsiPARP1                3.82e-07 ***
## poly(H2O2, 3)3:siRNAsiPARP1                0.906070    
## poly(H2O2, 3)1:genotypeALC1KO              0.079901 .  
## poly(H2O2, 3)2:genotypeALC1KO              3.71e-06 ***
## poly(H2O2, 3)3:genotypeALC1KO              0.362697    
## siRNAsiPARP1:genotypeALC1KO                8.14e-06 ***
## poly(H2O2, 3)1:siRNAsiPARP1:genotypeALC1KO 0.692050    
## poly(H2O2, 3)2:siRNAsiPARP1:genotypeALC1KO 1.61e-05 ***
## poly(H2O2, 3)3:siRNAsiPARP1:genotypeALC1KO 0.984250    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03987 on 32 degrees of freedom
## Multiple R-squared:  0.9922, Adjusted R-squared:  0.9886 
## F-statistic: 271.6 on 15 and 32 DF,  p-value: < 2.2e-16
cat("AIC: ", AIC(fit11))
## AIC:  -158.557
simres <- simulateResiduals(fittedModel = fit11)
plot(simres)

fit12 <- lmer(Counts ~ poly(H2O2,3)*siRNA*genotype + (1|UID), data = dataset)
print(summary(fit12))
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Counts ~ poly(H2O2, 3) * siRNA * genotype + (1 | UID)
##    Data: dataset
## 
## REML criterion at convergence: 439.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.67826 -0.45216 -0.01205  0.40378  1.92855 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  UID      (Intercept) 22071    148.6   
##  Residual             52429    229.0   
## Number of obs: 48, groups:  UID, 12
## 
## Fixed effects:
##                                             Estimate Std. Error        df
## (Intercept)                                 1067.833    108.288     8.000
## poly(H2O2, 3)1                             -4737.379    457.950    24.000
## poly(H2O2, 3)2                                61.116    457.950    24.000
## poly(H2O2, 3)3                              -743.243    457.950    24.000
## siRNAsiPARP1                                -230.583    153.142     8.000
## genotypeALC1KO                               -35.500    153.142     8.000
## poly(H2O2, 3)1:siRNAsiPARP1                 -544.032    647.639    24.000
## poly(H2O2, 3)2:siRNAsiPARP1                 1647.828    647.639    24.000
## poly(H2O2, 3)3:siRNAsiPARP1                    9.997    647.639    24.000
## poly(H2O2, 3)1:genotypeALC1KO              -1907.534    647.639    24.000
## poly(H2O2, 3)2:genotypeALC1KO               1754.587    647.639    24.000
## poly(H2O2, 3)3:genotypeALC1KO                106.013    647.639    24.000
## siRNAsiPARP1:genotypeALC1KO                  334.333    216.575     8.000
## poly(H2O2, 3)1:siRNAsiPARP1:genotypeALC1KO  -426.247    915.899    24.000
## poly(H2O2, 3)2:siRNAsiPARP1:genotypeALC1KO -1694.412    915.899    24.000
## poly(H2O2, 3)3:siRNAsiPARP1:genotypeALC1KO   -54.356    915.899    24.000
##                                            t value Pr(>|t|)    
## (Intercept)                                  9.861 9.42e-06 ***
## poly(H2O2, 3)1                             -10.345 2.53e-10 ***
## poly(H2O2, 3)2                               0.133  0.89495    
## poly(H2O2, 3)3                              -1.623  0.11766    
## siRNAsiPARP1                                -1.506  0.17057    
## genotypeALC1KO                              -0.232  0.82250    
## poly(H2O2, 3)1:siRNAsiPARP1                 -0.840  0.40919    
## poly(H2O2, 3)2:siRNAsiPARP1                  2.544  0.01780 *  
## poly(H2O2, 3)3:siRNAsiPARP1                  0.015  0.98781    
## poly(H2O2, 3)1:genotypeALC1KO               -2.945  0.00706 ** 
## poly(H2O2, 3)2:genotypeALC1KO                2.709  0.01225 *  
## poly(H2O2, 3)3:genotypeALC1KO                0.164  0.87134    
## siRNAsiPARP1:genotypeALC1KO                  1.544  0.16123    
## poly(H2O2, 3)1:siRNAsiPARP1:genotypeALC1KO  -0.465  0.64585    
## poly(H2O2, 3)2:siRNAsiPARP1:genotypeALC1KO  -1.850  0.07666 .  
## poly(H2O2, 3)3:siRNAsiPARP1:genotypeALC1KO  -0.059  0.95317    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("AIC: ", AIC(fit12))
## AIC:  475.7654
simres <- simulateResiduals(fittedModel = fit12)
plot(simres)

Compare Results

ICtab(fit1,fit2,fit3,fit4,
      fit5,fit6,fit7,fit8,
      fit9,fit10,fit11,fit12,
      base=T)
##       AIC    dAIC   df
## fit11 -158.6    0.0 17
## fit7  -119.1   39.5 13
## fit10  -73.8   84.7 17
## fit3   -65.6   92.9 9 
## fit6   -34.7  123.8 13
## fit2    24.3  182.9 9 
## fit12  475.8  634.3 18
## fit8   533.0  691.6 14
## fit4   621.3  779.8 10
## fit9   687.9  846.4 19
## fit5   689.8  848.4 15
## fit1   712.4  871.0 11

Final Result

fit <- fit11

output <- coef(summary(fit))
output <- output[grep("H2O2", rownames(output)),]


rownames(output) <- gsub("poly\\(|, [1-3]\\)","", rownames(output) )
rownames(output) <- gsub("genotype",  paste0(" ",levels(dataset$genotype)[1], " vs. "), rownames(output))
rownames(output)[!(grepl("vs", rownames(output)))] <- paste(rownames(output)[!(grepl("vs", rownames(output)))], levels(dataset$genotype)[1],  sep = " in " )

rownames(output) <- gsub("siRNA",  paste0(" ",levels(dataset$siRNA)[1], " vs. "), rownames(output))
rownames(output)[!(grepl("vs.*vs| in ", rownames(output)))] <- paste(rownames(output)[!(grepl("vs.*vs| in ", rownames(output)))], levels(dataset$siRNA)[1],  sep = " in " )

rownames(output)[!(grepl("vs", rownames(output)))] <- paste(rownames(output)[!(grepl("vs", rownames(output)))], levels(dataset$siRNA)[1],  sep = " " )


# suggested result table
kable(output, row.names = T)
Estimate Std. Error t value Pr(>|t|)
H2O21 in WT siCtrl -2.2700352 0.0797498 -28.4644590 0.0000000
H2O22 in WT siCtrl 0.0555822 0.0797498 0.6969573 0.4908621
H2O23 in WT siCtrl -0.3477220 0.0797498 -4.3601612 0.0001261
H2O21: siCtrl vs. siPARP1 in WT -0.1395743 0.1127833 -1.2375448 0.2248890
H2O22: siCtrl vs. siPARP1 in WT 0.7175920 0.1127833 6.3625752 0.0000004
H2O23: siCtrl vs. siPARP1 in WT 0.0134139 0.1127833 0.1189354 0.9060700
H2O21: WT vs. ALC1KO in siCtrl -0.2039924 0.1127833 -1.8087113 0.0799010
H2O22: WT vs. ALC1KO in siCtrl 0.6289899 0.1127833 5.5769793 0.0000037
H2O23: WT vs. ALC1KO in siCtrl 0.1041468 0.1127833 0.9234246 0.3626969
H2O21: siCtrl vs. siPARP1: WT vs. ALC1KO 0.0637482 0.1594996 0.3996763 0.6920503
H2O22: siCtrl vs. siPARP1: WT vs. ALC1KO -0.8089581 0.1594996 -5.0718495 0.0000161
H2O23: siCtrl vs. siPARP1: WT vs. ALC1KO 0.0031734 0.1594996 0.0198959 0.9842500
write.table(output, file = "FigureS3_H2O2_Stats_Ref_WT.txt", quote = F, sep = "\t", row.names = T, col.names = NA)
# re-fit with ALC1KO reference
dataset$genotype <- relevel(dataset$genotype, ref = "ALC1KO")


fit <- lm(NormCounts2 ~ poly(H2O2,3)*siRNA*genotype, data = dataset)

output <- coef(summary(fit))
output <- output[grep("H2O2", rownames(output)),]


rownames(output) <- gsub("poly\\(|, [1-3]\\)","", rownames(output) )
rownames(output) <- gsub("genotype",  paste0(" ",levels(dataset$genotype)[1], " vs. "), rownames(output))
rownames(output)[!(grepl("vs", rownames(output)))] <- paste(rownames(output)[!(grepl("vs", rownames(output)))], levels(dataset$genotype)[1],  sep = " in " )

rownames(output) <- gsub("siRNA",  paste0(" ",levels(dataset$siRNA)[1], " vs. "), rownames(output))
rownames(output)[!(grepl("vs.*vs| in ", rownames(output)))] <- paste(rownames(output)[!(grepl("vs.*vs| in ", rownames(output)))], levels(dataset$siRNA)[1],  sep = " in " )

rownames(output)[!(grepl("vs", rownames(output)))] <- paste(rownames(output)[!(grepl("vs", rownames(output)))], levels(dataset$siRNA)[1],  sep = " " )


# suggested result table
kable(output, row.names = T)
Estimate Std. Error t value Pr(>|t|)
H2O21 in ALC1KO siCtrl -2.4740275 0.0797498 -31.0223630 0.0000000
H2O22 in ALC1KO siCtrl 0.6845721 0.0797498 8.5839971 0.0000000
H2O23 in ALC1KO siCtrl -0.2435752 0.0797498 -3.0542416 0.0045206
H2O21: siCtrl vs. siPARP1 in ALC1KO -0.0758261 0.1127833 -0.6723172 0.5062076
H2O22: siCtrl vs. siPARP1 in ALC1KO -0.0913661 0.1127833 -0.8101032 0.4238645
H2O23: siCtrl vs. siPARP1 in ALC1KO 0.0165873 0.1127833 0.1470725 0.8839975
H2O21: ALC1KO vs. WT in siCtrl 0.2039924 0.1127833 1.8087113 0.0799010
H2O22: ALC1KO vs. WT in siCtrl -0.6289899 0.1127833 -5.5769793 0.0000037
H2O23: ALC1KO vs. WT in siCtrl -0.1041468 0.1127833 -0.9234246 0.3626969
H2O21: siCtrl vs. siPARP1: ALC1KO vs. WT -0.0637482 0.1594996 -0.3996763 0.6920503
H2O22: siCtrl vs. siPARP1: ALC1KO vs. WT 0.8089581 0.1594996 5.0718495 0.0000161
H2O23: siCtrl vs. siPARP1: ALC1KO vs. WT -0.0031734 0.1594996 -0.0198959 0.9842500
write.table(output, file = "FigureS3_H2O2_Stats_Ref_ALC1KO.txt", quote = F, sep = "\t", row.names = T, col.names = NA)

Anova

fit11a <- lm(NormCounts2 ~ poly(H2O2,3)*siRNA*genotype, data = dataset)
fit11b <- lm(NormCounts2 ~ poly(H2O2,3)*siRNA+genotype, data = dataset)

# anova table
anova(fit11a, fit11b)
## Analysis of Variance Table
## 
## Model 1: NormCounts2 ~ poly(H2O2, 3) * siRNA * genotype
## Model 2: NormCounts2 ~ poly(H2O2, 3) * siRNA + genotype
##   Res.Df     RSS Df Sum of Sq      F    Pr(>F)    
## 1     32 0.05088                                  
## 2     39 0.15961 -7  -0.10873 9.7692 1.861e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
fit11c <- lm(NormCounts2 ~ poly(H2O2,3)*genotype*siRNA, data = dataset)
fit11d <- lm(NormCounts2 ~ poly(H2O2,3)*genotype+siRNA, data = dataset)

# anova table
anova(fit11c, fit11d)
## Analysis of Variance Table
## 
## Model 1: NormCounts2 ~ poly(H2O2, 3) * genotype * siRNA
## Model 2: NormCounts2 ~ poly(H2O2, 3) * genotype + siRNA
##   Res.Df     RSS Df Sum of Sq      F   Pr(>F)    
## 1     32 0.05088                                 
## 2     39 0.16428 -7  -0.11339 10.188 1.21e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1